AccelerAge Bayesian Aft
Philippe Berends, 3972003
2025-06-12
1 Import libraries
2 Generation of survival data
generate_data = function(pop_n, obs_n, p, g, ltname, seed, force_recalc, X_plots,
X_rho, β_rho, non_zero_betas, X_scale, active_hazard,
betasep, bscale, betaplot = F, target_snr = 2,
method = "weibull", a, b) {
print("Generating βs...")
betas = generate_betas(p = p, g=g, rho = β_rho, rho_between = 0,
seed = seed, mu_u = betasep, mu_l = -betasep,
beta_scale = bscale, plot = betaplot,
non_zero_groups = non_zero_betas,
active_hazard = active_hazard,
target_snr = target_snr,
method = method, # Different scaling for different methods
weib_shape = a, weib_scale = b,
gomp_a = a, gomp_b = b)
# betas$beta = betas$beta + 0.01
print("Generating X...")
X = generate_X(n = pop_n, p = p, g = g, rho = X_rho, seed = seed,
rho_between = 0, X_plots = X_plots, scale = X_scale)
print("Done!")
if (method == "weibull") {
# a = shape
# b = scale
print("Generating Weibull population lifetable...")
lt = generate_population_lifetable_weibull(N_pop = pop_n, M = p,
betas = betas$beta, X = X$X,
filename = ltname, seed = seed,
force_recalc = force_recalc,
a = a, b = b) # Weibull params
print("Done!")
plot = ggplot(lt, aes(x = t, y = mrl)) +
geom_line() +
labs(title = "Weibull Population Lifetable MRL",
x = "Chronological age (years)", y = "Mean Residual Life") +
theme_minimal()
print(plot)
print("Generating data from Weibull lifetable...")
df_sim = create_dataset_weibull(M = p, n_obs = obs_n, G = g,
gsize = (p/g), lt = lt, betas = betas$beta,
seednr = seed+1, X_plots = F,
a = a, b = b, # Weibull params
X_rho = X_rho, X_scale = X_scale,
followup = 20) # Nr years followed
} else {
print("Generating Gompertz population lifetable...")
lt = generate_population_lifetable_gompertz(N_pop = pop_n, M = p,
betas = betas$beta, X = X$X,
filename = ltname, seed = seed,
force_recalc = force_recalc,
a = a, b = b) # Gompertz params
print("Done!")
plot = ggplot(lt, aes(x = t, y = mrl)) +
geom_line() +
labs(title = "Gompertz Population Lifetable MRL",
x = "Chronological age (years)", y = "Mean Residual Life") +
theme_minimal()
print(plot)
print("Generating data from Gompertz lifetable...")
df_sim = create_dataset_gompretz(M = p, n_obs = obs_n, G = g,
gsize = (p/g), lt = lt, betas = betas$beta,
seednr = seed+1, X_plots = F,
a = a, b = b, # Gompertz params
X_rho = X_rho, X_scale = X_scale,
followup = 20) # Nr years followed
}
print("Done!")
return(list(df = df_sim, true_betas = betas, lt = lt))
}3 Functions to estimate MRL
3.1 Frequentist AFT
aft_reg = function(df_sim, p, method) {
result = switch(method,
"gompertz" = aft_reg_gompertz(df_sim, p),
"weibull" = aft_reg_weibull(df_sim, p),
stop("Unsupported method: ", method)
)
return(result)
}3.1.1 Gompertz implementation
aft_reg_gompertz = function(df_sim, p) {
true_betas = df_sim$true_betas
lt = df_sim$lt
df_sim = df_sim$df
# Train/test split 50/50
train_indices = sample(1:nrow(df_sim), nrow(df_sim)/2)
df_sim_train = df_sim[train_indices,]
df_sim_test = df_sim[-train_indices,]
# Does not work in highdim cases
if (dim(df_sim_train)[1] < p) {
return(NULL)
}
#### FIT MODEL
predictor_vars = paste(c(glue("x{1:p}")), collapse = " + ")
aft_formula = as.formula(paste("Surv(age_start, age_end, status) ~", predictor_vars))
fit_aft = tryCatch({
eha::aftreg(
formula = aft_formula,
data = df_sim_train,
dist = "gompertz",
param = "lifeExp",
control = list(maxit = 100, trace = FALSE)
)
}, error = function(e) {
message("AFT model error: ", e$message)
NULL
})
if (is.null(fit_aft)) {
print("Fit failed")
return(NULL)
}
sigma_est <- exp(fit_aft$coefficients["log(scale)"])
tau_est <- exp(fit_aft$coefficients["log(shape)"])
a_est <- tau_est / sigma_est
b_est <- 1 / sigma_est
betas_est_aft <- fit_aft$coefficients[1:p]
linpred_aft <- rowSums(sweep(df_sim_test[,1:p], 2, betas_est_aft, "*"))
distr_ests = list(sigma = sigma_est, tau = tau_est, a = a_est, b = b_est)
#### ESTIMATE MRL
df_sim_test$pred_mrl = NA
for (i in 1:nrow(df_sim_test)) {
mrl_uncon <- integrate(gomp_baseline_surv, lower = (df_sim_test$age_start[i] * exp(linpred_aft[i])),
upper=Inf, a = a_est, b = b_est)$value
s_cond <- gomp_baseline_surv(df_sim_test$age_start[i] * exp(linpred_aft[i]), a = a_est, b = b_est)
df_sim_test$pred_mrl[i] <- (mrl_uncon / s_cond) * exp(-linpred_aft[i])
if(is.na(df_sim_test$pred_mrl[i]) || df_sim_test$pred_mrl[i] < 0) {
df_sim_test$pred_mrl[i] <- NA
}
}
pred_mrl = df_sim_test$pred_mrl
# obtain biological age, assuming true lifetable is known
aft_b = vector(length = length(pred_mrl))
for (i in 1:length(pred_mrl)){
aft_b[i] = lt$t[which.min(abs(lt$mrl - pred_mrl[i]))]
}
# Calculate rmses
aftreg_coef_rmse = sqrt(mean((betas_est_aft - true_betas$beta)^2))
aftreg_rmse = sqrt(mean((df_sim_test$mrl - df_sim_test$pred_mrl)^2))
aftreg_bioage_rmse = sqrt(mean((aft_b - df_sim_test$b_age)^2))
# plt = df_sim_test %>% ggplot(aes(x = mrl, y = pred_mrl)) +
# geom_point() +
# geom_abline(color = 'red') +
# labs(title = "True vs Predicted MRL for EHA AFT Gompertz",
# subtitle = glue("RMSE = {aftreg_rmse}")) +
# theme_minimal()
# print(plt)
true_betas$betahat = betas_est_aft
# plt = true_betas %>%
# mutate(index = 1:p) %>%
# ggplot(aes(x = index)) +
# geom_segment(aes(xend = index, y = beta, yend = betahat), color = "gray") +
# geom_point(aes(y = beta, color = "True")) +
# geom_point(aes(y = betahat, color = "Predicted")) +
# scale_color_manual(values = c("True" = "black", "Predicted" = "skyblue"),
# name = "Type") +
# ylab("β") +
# xlab("Index") +
# labs(title = "True vs Predicted β-coefficients for AFT Gompertz",
# subtitle = glue("RMSE = {aftreg_coef_rmse %>% round(4)}")) +
# theme_minimal()
# print(plt)
return(list(
test_data = df_sim_test,
true_betas = true_betas,
predicted_betas = betas_est_aft,
predicted_bioage = aft_b,
rmses = list(coef = aftreg_coef_rmse, mrl = aftreg_rmse, bioage = aftreg_bioage_rmse),
method = "AFT gompertz",
distr_ests = distr_ests
))
}3.1.2 Weibull implementation
aft_reg_weibull = function(df_sim, p) {
true_betas = df_sim$true_betas
lt = df_sim$lt
df_sim = df_sim$df
# Train/test split 50/50
train_indices = sample(1:nrow(df_sim), nrow(df_sim)/2)
df_sim_train = df_sim[train_indices,]
df_sim_test = df_sim[-train_indices,]
# Does not work in highdim cases
if (dim(df_sim_train)[1] < p) {
return(NULL)
}
#### FIT MODEL
predictor_vars = paste(c(glue("x{1:p}")), collapse = " + ")
aft_formula = as.formula(paste("Surv(age_start, age_end, status) ~", predictor_vars))
fit_aft = tryCatch({
eha::aftreg(
formula = aft_formula,
data = df_sim_train,
dist = "weibull",
param = "lifeExp", # Otherwise beta signs are inverted
control = list(maxit = 100, trace = FALSE)
)
}, error = function(e) {
message("AFT model error: ", e$message)
NULL
})
fa <<- fit_aft
if (is.null(fit_aft)) {
print("Fit failed")
return(NULL)
}
print("----------")
shape_est <- exp(fit_aft$coefficients["log(shape)"])
cat("shape_est:", shape_est, "\n")
scale_est <- exp(fit_aft$coefficients["log(scale)"])
cat("scale_est:", scale_est, "\n")
betas_est_aft <- fit_aft$coefficients[1:p]
linpred_aft <- rowSums(sweep(df_sim_test[,1:p], 2, betas_est_aft, "*"))
# scale^(-shape) == lambda so we use it as est (from Marije's script)
lambda_est <- scale_est^(-shape_est)
nu_est <- shape_est
#### ESTIMATE MRL
df_sim_test$pred_mrl = NA
for (i in 1:nrow(df_sim_test)) {
scale_adj = scale_est * exp(linpred_aft[i])
mrl_uncon <- integrate(function(t) {
1 - pweibull(t, shape = shape_est, scale = scale_adj)
}, lower = df_sim_test$age_start[i], upper = Inf)$value
s_cond <- 1 - pweibull(df_sim_test$age_start[i], shape = shape_est, scale = scale_adj)
df_sim_test$pred_mrl[i] <- mrl_uncon / s_cond
# mrl_uncon <- integrate(weib_baseline_surv,
# lower = (df_sim_test$age_start[i] * exp(linpred_aft[i])),
# upper = Inf,
# lambda = lambda_est, nu = nu_est)$value
# s_cond <- weib_baseline_surv(df_sim_test$age_start[i] * exp(linpred_aft[i]),
# lambda = lambda_est, nu = nu_est)
# df_sim_test$pred_mrl[i] <- (mrl_uncon / s_cond) * exp(-linpred_aft[i])
}
print("----------")
# obtain biological age, assuming true lifetable is known
pred_mrl = df_sim_test$pred_mrl
aft_b = vector(length = length(pred_mrl))
for (i in 1:length(pred_mrl)){
aft_b[i] = lt$t[which.min(abs(lt$mrl - pred_mrl[i]))]
}
# Calculate rmses
aftreg_coef_rmse = sqrt(mean((betas_est_aft - true_betas$beta)^2))
aftreg_rmse = sqrt(mean((df_sim_test$mrl - df_sim_test$pred_mrl)^2))
aftreg_bioage_rmse = sqrt(mean((aft_b - df_sim_test$b_age)^2))
plt = df_sim_test %>% ggplot(aes(x = mrl, y = pred_mrl)) +
geom_point() +
geom_abline(color = 'red') +
labs(title = "True vs Predicted MRL for EHA AFT Weibull",
subtitle = glue("RMSE = {aftreg_rmse}")) +
theme_minimal()
# print(plt)
true_betas$betahat = betas_est_aft
plt = true_betas %>%
mutate(index = 1:p) %>%
ggplot(aes(x = index)) +
geom_segment(aes(xend = index, y = beta, yend = betahat), color = "gray") +
geom_point(aes(y = beta, color = "True")) +
geom_point(aes(y = betahat, color = "Predicted")) +
scale_color_manual(values = c("True" = "black", "Predicted" = "skyblue"),
name = "Type") +
ylab("β") +
xlab("Index") +
labs(title = "True vs Predicted β-coefficients for AFT Weibull",
subtitle = glue("RMSE = {aftreg_coef_rmse %>% round(4)}")) +
theme_minimal()
# print(plt)
return(list(
test_data = df_sim_test,
true_betas = true_betas,
predicted_betas = betas_est_aft,
predicted_bioage = aft_b,
rmses = list(coef = aftreg_coef_rmse, mrl = aftreg_rmse, bioage = aftreg_bioage_rmse),
method = "AFT weibull",
distr_ests = list(shape = shape_est, scale = scale_est)
))
}3.2 PSBC Bayesian AFT
psbc_reg = function(df_sim, p, g) {
true_betas = df_sim$true_betas
df_sim = df_sim$df
lt = df_sim$lt
# Train/test split 50/50
train_indices = sample(1:nrow(df_sim), nrow(df_sim)/2)
df_sim_train = df_sim[train_indices,]
df_sim_test = df_sim[-train_indices,]
Y_train = cbind(df_sim_train$age_start, df_sim_train$age_end, df_sim_train$status)
Y_test = cbind(df_sim_test$age_start, df_sim_test$age_end, df_sim_test$status)
# Create covariate matrix X
X_train = as.matrix(df_sim_train[, paste0("x", 1:p)])
X_test = as.matrix(df_sim_test[, paste0("x", 1:p)])
XC = NULL # CONFOUNDERS
# Create group indicator for each variable
grpInx = rep(1:g, each = p/g)
# Set hyperparameters (PRIOR)?
hyperParams = list(
a.sigSq = 2, b.sigSq = 2,
mu0 = 0, h0 = 10^6,
v = 10^6
)
# Set starting values https://cran.r-project.org/web/packages/psbcGroup/psbcGroup.pdf
startValues = list(
w = log(Y_train[,1]),
beta = runif(p, -0.5, 0.5), # uniform draws from realistic range
tauSq = rep(0.1, p),
mu = 0.1,
sigSq = 0.1,
lambdaSq = 0.1,
betaC = rep(0.11, 0)
)
# MCMC parameters format __ CHANGE __
mcmcParams = list(
run = list(
numReps = 5000,
thin = 5, # Update params every n steps
burninPerc = 0.2 # Discard the first ... steps
),
tuning = list(
mu.prop.var = 0.1,
sigSq.prop.var = 0.005,
L.beC = 50,
M.beC = 1,
eps.beC = 0.001,
L.be = 100,
M.be = 1,
eps.be = 0.001
)
)
fit_bayes_aft = aftGL_LT(Y_train, X_train, XC, grpInx, hyperParams, startValues, mcmcParams)
beta_samples = fit_bayes_aft$beta.p
# # Convert to long format
# beta_long = beta_samples %>%
# as.data.frame() %>%
# mutate(iteration = row_number()) %>%
# pivot_longer(cols = -iteration,
# names_to = "beta_index",
# values_to = "value") %>%
# mutate(beta_index = factor(beta_index, levels = paste0("V", 1:ncol(beta_samples))))
# # Create boxplot
# plt = ggplot(beta_long, aes(x = beta_index, y = value)) +
# geom_boxplot(fill = "lightblue", alpha = 0.7) +
# geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
# labs(
# title = "Posterior Distributions of Beta Coefficients",
# x = "Beta Index",
# y = "Coefficient Value"
# ) +
# theme_minimal() +
# theme(axis.text.x = element_text(angle = 45, hjust = 1))
#
# print(plt)
# Extract betas using the maximum likelihood via density
get_mode = function(x){
xdist = density(x)
mode = xdist$x[which.max(xdist$y)]
return(mode)
}
# betas = fit_bayes_aft$beta.p %>% apply(MARGIN = 2, FUN = get_mode)
betas = apply(fit_bayes_aft$beta.p, 2, mean)
# Make linear predictor for MRL
X_test = X_test %>% as.matrix
linpred = X_test %*% betas
# Change to lognormal
# Create baseline lognormal survival
# Implemented in R -> check log scales
# dlnorm, use any of them 1 - F(x)
# instead of gomp_baseline_surv no false, 1-
# a = exp(-9)
# b = 0.085
# 1-plnorm(1:100, meanlog = fit_bayes_aft$mu.p[120], sdlog = fit_bayes_aft$sigSq.p[120])
# Estimate mean from MCMC draws
mu_est = mean(fit_bayes_aft$mu.p)
sigSq_est = mean(fit_bayes_aft$sigSq.p)
for (i in 1:nrow(df_sim_test)){
mrl_uncon = integrate(function(x) 1 - plnorm(x, meanlog = mu_est, sdlog = sqrt(sigSq_est)), # Var to Sd via sqrt
lower = (df_sim_test$age_start[i] * exp(linpred[i])),
upper = Inf)$value
s_cond = 1 - plnorm(df_sim_test$age_start[i] * exp(linpred[i]), meanlog = mu_est, sdlog = sqrt(sigSq_est))
df_sim_test$pred_mrl[i] <- (mrl_uncon / s_cond) * exp(-linpred[i])
}
# Get the RMSE
pred_RMSE = sqrt(mean((df_sim_test$mrl - df_sim_test$pred_mrl)^2))
beta_RMSE =sqrt(mean((betas - true_betas$beta)^2))
# MRL prediction plot
plt = df_sim_test %>% ggplot(aes(x = mrl, y = pred_mrl)) +
geom_point() +
geom_abline(color = 'red') +
labs(title = "True vs Predicted MRL for PSBC Bayesian AFT",
subtitle = glue("RMSE = {pred_RMSE}")) +
theme_minimal()
# print(plt)
# Plot difference between true and estimates betas
true_betas$betahat = betas
plt = true_betas %>%
mutate(index = 1:p) %>%
ggplot(aes(x = index)) +
# Vertical lines between pred and actual
geom_segment(aes(xend = index, y = beta, yend = betahat), color = "gray") +
geom_point(aes(y = beta, color = "True")) +
geom_point(aes(y = betahat, color = "Predicted")) +
scale_color_manual(values = c("True" = "black", "Predicted" = "skyblue"),
name = "Type") +
ylab("β") +
xlab("Index") +
labs(title = "True vs Predicted β-coefficients for Bayesian AFT", subtitle = glue("RMSE = {beta_RMSE %>% round(4)}")) +
theme_minimal()
# print(plt)
return(list(
test_data = df_sim_test,
true_betas = true_betas,
predicted_betas = betas,
method = "PSBC",
distr_ests = list(mu = mu_est, sigSq = sigSq_est)
))
}3.3 BRMS Bayesian AFT
brms_reg = function(df_sim, p, g, method, seed = 123) {
result = switch(method,
"gompertz" = brms_reg_gompertz(df_sim, p, g, seed),
"weibull" = brms_reg_weibull(df_sim, p, g, seed),
stop("Unsupported method: ", method)
)
return(result)
}3.3.1 Gompertz implementation
brms_reg_gompertz = function(df_sim, p, g, seed) {
# Grompertz stan implementation source:
# https://metrumrg.com/wp-content/uploads/2024/12/2024-Yoder-cure-rate-brms.pdf
# Todd Yoder, Andrew Tredennick, Timothy Waterhouse: Gompertz Cure Rate Survival Models with Stan and brms
stan_funs <- "
real gompertz_lpdf(real t) {
return log(mu) + gamma*t + exp(gamma*t);
}
real gompertz_lcdf(real t) {
return log(1 - exp((1 - exp(gamma*t))));
}
real gompertz_lccdf(real t) {
return (1 - exp(gamma*t));
}
real gompertz_rng(real mu, real gamma) {
real sim_time;
sim_time = 1/gamma*log(1 - gamma/mu*log(1 - uniform_rng(0, 1)));
if (is_nan(sim_time))
sim_time = positive_infinity();
return sim_time;
}
"
gompertz_fam <- brms::custom_family(
name = "gompertz",
dpars = c("mu", "gamma"),
links = c("log", "identity"),
lb = c(0, NA),
type = "real",
log_lik = function(i, prep) {
mu <- brms::get_dpar(prep, "mu", i = i)
gamma <- brms::get_dpar(prep, "gamma", i = i)
t <- prep$data$Y[i]
cens <- prep$data$cens[i]
if (cens == 0) x <- gompertz_lpdf(t, mu, gamma)
if (cens == 1) x <- gompertz_lccdf(t, mu, gamma)
return(x)
},
posterior_predict = function(i, prep, ...) {
mu <- brms::get_dpar(prep, "mu", i = i)
gamma <- brms::get_dpar(prep, "gamma", i = i)
return(gompertz_rng(mu, gamma))
}
)
true_betas = df_sim$true_betas
lt = df_sim$lt
df_sim = df_sim$df
# Train/test split 50/50
train_indices = sample(1:nrow(df_sim), nrow(df_sim)/2)
df_sim_train = df_sim[train_indices,]
df_sim_test = df_sim[-train_indices,]
#### FIT MODEL
preds = paste(paste0("x", 1:p), collapse = " + ")
formula = as.formula(paste("age_end | cens(1 - status) + trunc(lb = age_start) ~", preds))
priors = c(
prior(normal(0, 0.5), class = "b"), # BETAS, using ridge, λ = 2 in this case
prior(normal(log(80), 0.4), class = Intercept)#, # about human life expectancy
# prior(gamma(2,1), class = shape) # Hazard
)
fit_bayes_aft = brm(
formula = formula,
data = df_sim_train,
family = gompertz_fam,
stanvars = stanvar(scode = stan_funs, block = "functions"),
prior = priors,
warmup = 2000, iter = 4000, chains = 4, cores = 4
)
# Extract posterior samples and calculate estimates
post = posterior_samples(fit_bayes_aft)
beta_cols = paste0("b_x", 1:p)
betahat = apply(post[, beta_cols], 2, mean)
X_test = df_sim_test[,1:p] %>% as.matrix
linpred = X_test %*% betahat
# Extract Gompertz parameters from posterior
mu_est = mean(post$mu)
gamma_est = mean(post$gamma)
# Convert to a/b parameterization for consistency with AFT
a_est = gamma_est
b_est = 1/mu_est
#### ESTIMATE MRL (using median approach like AFT)
df_sim_test$pred_mrl = NA
for (i in 1:nrow(df_sim_test)) {
mrl_uncon <- integrate(gomp_baseline_surv, lower = (df_sim_test$age_start[i] * exp(linpred[i])),
upper=Inf, a = a_est, b = b_est)$value
s_cond <- gomp_baseline_surv(df_sim_test$age_start[i] * exp(linpred[i]), a = a_est, b = b_est)
df_sim_test$pred_mrl[i] <- (mrl_uncon / s_cond) * exp(-linpred[i])
if(is.na(df_sim_test$pred_mrl[i]) || df_sim_test$pred_mrl[i] < 0) {
df_sim_test$pred_mrl[i] <- NA
}
}
# obtain biological age, assuming true lifetable is known
pred_mrl = df_sim_test$pred_mrl
brms_b = vector(length = length(pred_mrl))
for (i in 1:length(pred_mrl)){
brms_b[i] = lt$t[which.min(abs(lt$mrl - pred_mrl[i]))]
}
# Calculate rmses
pred_RMSE = sqrt(mean((df_sim_test$mrl - df_sim_test$pred_mrl)^2))
beta_RMSE = sqrt(mean((betahat - true_betas$beta)^2))
bioage_RMSE = sqrt(mean((brms_b - df_sim_test$b_age)^2))
# MRL prediction plot
plt = df_sim_test %>% ggplot(aes(x = mrl, y = pred_mrl)) +
geom_point() +
geom_abline(color = 'red') +
labs(title = "True vs Predicted MRL for BRMS AFT Gompertz",
subtitle = glue("RMSE = {pred_RMSE}")) +
theme_minimal()
# print(plt)
# Beta comparison plot
true_betas$betahat = betahat
plt = true_betas %>%
mutate(index = 1:p) %>%
ggplot(aes(x = index)) +
geom_segment(aes(xend = index, y = beta, yend = betahat), color = "gray") +
geom_point(aes(y = beta, color = "True")) +
geom_point(aes(y = betahat, color = "Predicted")) +
scale_color_manual(values = c("True" = "black", "Predicted" = "skyblue"),
name = "Type") +
ylab("β") +
xlab("Index") +
labs(title = "True vs Predicted β-coefficients for BRMS AFT Gompertz",
subtitle = glue("RMSE = {beta_RMSE %>% round(4)}")) +
theme_minimal()
# print(plt)
return(list(
test_data = df_sim_test,
true_betas = true_betas,
predicted_betas = betahat,
predicted_bioage = brms_b,
rmses = list(coef = beta_RMSE, mrl = pred_RMSE, bioage = bioage_RMSE),
method = "BRMS gompertz",
distr_ests = list(mu = mu_est, gamma = gamma_est, a = a_est, b = b_est)
))
}3.3.2 Weibull implementation
brms_reg_weibull = function(df_sim, p, g, seed = 123) {
true_betas = df_sim$true_betas
lt = df_sim$lt
df_sim = df_sim$df
# Train/test split 50/50
train_indices = sample(1:nrow(df_sim), nrow(df_sim)/2)
df_sim_train = df_sim[train_indices,]
df_sim_test = df_sim[-train_indices,]
#### FIT MODEL
preds = paste(paste0("x", 1:p), collapse = " + ")
formula = as.formula(paste("age_end | cens(1 - status) + trunc(lb = age_start) ~", preds))
priors = c(
#prior(normal(0,0.5), class = "b"), # I took a wide sd due to convergence problems
prior(horseshoe(df = 3, par_ratio = 0.1), class = "b"), # WOW!!!
prior(normal(log(80), 0.2), class = Intercept), # about human life expectancy
# prior(normal(10, 4), class = shape) # Hazard: bad, gamma better?
prior(gamma(50, 5), class = shape) # not log scale!! (mean of gamma is a/b)
)
fit_bayes_aft = brm(
formula = formula,
data = df_sim_train,
family = weibull(),
prior = priors,
# warmup = 250,
iter = 6000, chains = 2, cores = 2,
control = list(adapt_delta = 0.97, max_treedepth = 15, stepsize = 0.01),
seed = seed,
# install.packages("cmdstanr", repos = c('https://stan-dev.r-universe.dev', getOption("repos")))
backend = "cmdstanr",
threads = threading(2),
save_pars = save_pars(all = TRUE), init = 0
)
# Whenever you see the warning
# "There were x divergent transitions after warmup.", you should really think about
# increasing adapt_delta. To do this, write control = list(adapt_delta = <x>), where
# <x> should usually be a value between 0.8 (current default) and 1. Increasing adapt_delta
# will slow down the sampler but will decrease the number of divergent transitions threatening
# the validity of your posterior samples.
# https://cran.r-project.org/web/packages/brms/vignettes/brms_overview.pdf
fb <<- fit_bayes_aft
# Extract posterior samples and calculate estimates
post = posterior_samples(fit_bayes_aft)
beta_cols = paste0("b_x", 1:p)
betahat = apply(post[, beta_cols], 2, mean)
X_test = df_sim_test[,1:p] %>% as.matrix
linpred = X_test %*% betahat
shape_est = (mean(post$shape)) # Not exp because not in log scale????
cat("Shape estimate:", shape_est, "\n")
scale_est = exp(mean(post$b_Intercept))
cat("Scale estimate", (scale_est), "\n")
#### ESTIMATE MRL
# scale^(-shape) == lambda so we use it as est (from marijes script)
lambda_est <- scale_est^(-shape_est)
nu_est <- shape_est
for (i in 1:nrow(df_sim_test)) {
scale_adj = scale_est * exp(linpred[i])
mrl_uncon <- integrate(function(t) {
1 - pweibull(t, shape = shape_est, scale = scale_adj)
}, lower = df_sim_test$age_start[i], upper = Inf)$value
# Conditional survival at age_start
s_cond <- 1 - pweibull(df_sim_test$age_start[i], shape = shape_est, scale = scale_adj)
# Mean residual life
df_sim_test$pred_mrl[i] <- mrl_uncon / s_cond
# mrl_uncon <- integrate(weib_baseline_surv,
# lower = (df_sim_test$age_start[i] * exp(linpred[i])),
# upper=Inf, lambda = lambda_est, nu = nu_est)$value
# s_cond <- weib_baseline_surv(df_sim_test$age_start[i] * exp(linpred[i]),
# lambda = lambda_est, nu = nu_est)
# df_sim_test$pred_mrl[i] <- (mrl_uncon / s_cond) * exp(-linpred[i])
}
# obtain biological age, assuming true lifetable is known
pred_mrl = df_sim_test$pred_mrl
brms_b = vector(length = length(pred_mrl))
for (i in 1:length(pred_mrl)){
brms_b[i] = lt$t[which.min(abs(lt$mrl - pred_mrl[i]))]
}
# Calculate rmses
pred_RMSE = sqrt(mean((df_sim_test$mrl - df_sim_test$pred_mrl)^2))
beta_RMSE = sqrt(mean((betahat - true_betas$beta)^2))
bioage_RMSE = sqrt(mean((brms_b - df_sim_test$b_age)^2))
# MRL prediction plot
# plt = df_sim_test %>% ggplot(aes(x = mrl, y = pred_mrl)) +
# geom_point() +
# geom_abline(color = 'red') +
# labs(title = "True vs Predicted MRL for BRMS AFT Weibull",
# subtitle = glue("RMSE = {pred_RMSE}")) +
# theme_minimal()
# print(plt)
# Beta comparison plot
true_betas$betahat = betahat
# plt = true_betas %>%
# mutate(index = 1:p) %>%
# ggplot(aes(x = index)) +
# geom_segment(aes(xend = index, y = beta, yend = betahat), color = "gray") +
# geom_point(aes(y = beta, color = "True")) +
# geom_point(aes(y = betahat, color = "Predicted")) +
# scale_color_manual(values = c("True" = "black", "Predicted" = "skyblue"),
# name = "Type") +
# ylab("β") +
# xlab("Index") +
# labs(title = "True vs Predicted β-coefficients for BRMS AFT Weibull",
# subtitle = glue("RMSE = {beta_RMSE %>% round(4)}")) +
# theme_minimal()
# print(plt)
return(list(
test_data = df_sim_test,
true_betas = true_betas,
predicted_betas = betahat,
predicted_bioage = brms_b,
rmses = list(coef = beta_RMSE, mrl = pred_RMSE, bioage = bioage_RMSE),
method = "BRMS weibull",
distr_ests = list(shape = shape_est, scale = scale_est)
))
}3.4 Monte Carlo Experiments
run_grid_experiments_mc = function(
n_grid = c(200, 500, 1000),
p_grid = c(20, 60, 100),
g_grid = c(4, 6, 10),
non_zero_grid = c(0.5, 0.75, 1.0),
M = 10,
pop_n = 100000,
bsep = 4,
X_rho = 0,
target_snr = 0.1,
method = "weibull",
seed_base = 123
) {
# Create parameter grid
param_grid = expand.grid(
n_obs = n_grid,
p = p_grid,
g = g_grid,
non_zero = non_zero_grid
)
# Filter invalid combinations, p should be divisible by g (not irrational)
param_grid = param_grid[param_grid$p %% param_grid$g == 0, ]
cat(glue("Running {nrow(param_grid)} parameter combinations with {M} replications each\n"))
cat(glue("Total experiments: {nrow(param_grid) * M}\n\n"))
all_results = list()
if (method == "gompertz") {
a = exp(-9)
b = 0.085
} else if (method == "weibull") {
a = 5 # Shape
b = 80 # Scale
}
# Iterate through the grid
for (i in 1:nrow(param_grid)) {
# Params to compare
params = param_grid[i,]
cat(glue("\n<<<<<<<< Parameter Set {i}/{nrow(param_grid)} >>>>>>>>>>"), "\n")
cat(glue("n = {params$n_obs} | p = {params$p} | g = {params$g} | non-zero = {params$non_zero}"), "\n")
# Run M Monte Carlo replications
mc_results = list()
for (m in 1:M) {
seed = seed_base + m
cat(glue(" Rep {m}/{M}... "), "\n")
# Generate data
gres = generate_data(
pop_n = pop_n,
obs_n = params$n_obs,
p = params$p,
g = params$g,
ltname = glue("mc_{.....}", "\n"),
seed = seed,
force_recalc = FALSE, # IMPORTANT! SAME LT!
X_plots = FALSE,
X_rho = X_rho,
β_rho = 0.9,
non_zero_betas = params$non_zero,
X_scale = 1,
active_hazard = 0,
betasep = bsep,
bscale = 1,
betaplot = FALSE,
target_snr = target_snr,
method = method,
a = a,
b = b
)
results = list()
# AFT
print("Running Frequentist AFT...")
set.seed(seed)
results$aft = tryCatch(
aft_reg(df_sim = gres, p = params$p, method = method),
error = function(e) NULL
)
# # PSBC
# print("Running PSBC Bayesian AFT...")
# set.seed(seed)
# results$psbc = tryCatch(
# psbc_reg(df_sim = df_sim, p = params$p, g = params$g),
# error = function(e) NULL
# )
# BRMS
print("Running BRMS Bayesian AFT...")
set.seed(seed)
results$brms = tryCatch(
brms_reg(df_sim = gres, p = params$p, g = params$g, method = method),
error = function(e) NULL
)
mc_results[[m]] = list(
params = params,
seed = seed,
results = results
)
cat("Done\n")
}
## TODO: SAVE RESULTS TO FILES IN SOME EXPERIMENT FOLDER
all_results[[i]] = list(
params = params,
mc_results = mc_results
)
}
return(all_results)
}4 Experiments Monte Carlo
4.1 \(n \in \{200, 400, 1000\}, p \in \{20\}, g = 5\)
5 High dimensional example (p = n)
Since obs_n = 200, train test split of 50/50 means that it is trained on 100 observations.
# rmses = run_experiment(p = 10, g = 4, ltname = "p60_may20", seed = 43, nzb = 0.75, bsep = 4, X_rho = 0, obs_n = 200)
# print(rmses)Small tests:
6 DEBUG - before performing actual MC experiments
6.1 Debug function:
debug_methods = function(n = 100, p = 10, g = 5, method = "weibull", seed = 123, t_SNR = 2, perf_A = T, nzg = 1) {
cat("**** SETTINGS ****\n")
cat("n =", n, "p =", p, "g =", g, "method =", method, "\n\n")
cat("Generating data...\n")
if (method == "weibull") {
a = 10 # shape
b = 80 # scale
} else {
a = exp(-9) # gompertz a
b = 0.085 # gompertz b
}
gres = generate_data(
pop_n = as.integer(1e4), # 10,000
obs_n = n * 2, # x2 due to train test split
p = p,
g = g,
ltname = glue("debug_{method}"),
seed = seed,
force_recalc = TRUE,
X_plots = FALSE,
X_rho = 0,
β_rho = 0.9,
non_zero_betas = nzg,
X_scale = 1,
active_hazard = 0,
betasep = 4,
bscale = 1,
betaplot = TRUE,
target_snr = t_SNR,
method = method,
a = a,
b = b
)
cat("Data generated successfully!\n")
cat("Dataset size:", nrow(gres$df), "x", ncol(gres$df), "\n\n")
if (!perf_A) {
return(gres)
}
cat("**** TESTING AFT ****\n")
set.seed(seed)
aft_result = aft_reg(df_sim = gres, p = p, method = method)
if (!is.null(aft_result)) {
cat("AFT SUCCESS!\n")
cat("AFT RMSEs:", "\n")
print(aft_result$rmses)
plot(fa)
} else {
cat("AFT FAILED!\n")
}
# aft_result = NULL # DEBUG
cat("\n**** TESTING BRMS ****\n")
set.seed(seed)
brms_result = brms_reg(df_sim = gres, p = p, g = g, method = method, seed = seed)
if (!is.null(brms_result)) {
cat("BRMS SUCCESS!\n")
cat("BRMS RMSEs:", "\n")
print(brms_result$rmses)
plot(fb, variable = c("b_Intercept", "shape", "b_x1"))
} else {
cat("BRMS FAILED!\n")
}
res = list(aft = aft_result, brms = brms_result)
cat("\n**** COMPARISON ****\n")
if (!is.null(brms_result)) {
if (is.null(aft_result)) {
res$aft = list(
test_data = data.frame(
mrl = rep(0, n),
pred_mrl = rep(0, n),
b_age = rep(0, n)
),
true_betas = res$brms$true_betas$beta,
predicted_betas = rep(0, p),
predicted_bioage = rep(0, n),
rmses = list(coef = 0, mrl = 0, bioage = 0),
method = "AFT (failed)"
)
aft_result = res$aft
}
p_vars = length(res$brms$true_betas$beta)
beta_comparison = data.frame(
index = 1:p_vars,
true_beta = res$brms$true_betas$beta,
aft_beta = res$aft$predicted_betas,
brms_beta = res$brms$predicted_betas
)
# Beta plot
p1 = ggplot(beta_comparison, aes(x = index)) +
geom_segment(aes(xend = index, y = true_beta, yend = aft_beta),
color = "lightgray", alpha = 0.5) +
geom_segment(aes(xend = index, y = true_beta, yend = brms_beta),
color = "lightblue", alpha = 0.5) +
geom_point(aes(y = true_beta, color = "True"), size = 1.5) +
geom_point(aes(y = aft_beta, color = "AFT"), size = 1) +
geom_point(aes(y = brms_beta, color = "BRMS"), size = 1) +
scale_color_manual(values = c("True" = "black", "AFT" = "red", "BRMS" = "blue")) +
geom_hline(yintercept = 0) +
labs(title = "Beta Coefficient Comparison: AFT vs BRMS",
subtitle = paste("AFT Beta RMSE:", round(res$aft$rmses$coef, 4),
"| BRMS Beta RMSE:", round(res$brms$rmses$coef, 4)),
x = "Coefficient Index", y = "Beta Value") +
theme_classic()
print(p1)
n_test = min(nrow(res$aft$test_data), nrow(res$brms$test_data))
mrl_comparison = data.frame(
true_mrl = res$brms$test_data$mrl[1:n_test],
aft_mrl = res$aft$test_data$pred_mrl[1:n_test],
brms_mrl = res$brms$test_data$pred_mrl[1:n_test]
)
mrl_comparison = mrl_comparison[complete.cases(mrl_comparison), ]
p2 = ggplot(mrl_comparison, aes(x = true_mrl)) +
geom_point(aes(y = aft_mrl, color = "AFT"), alpha = 0.6) +
geom_point(aes(y = brms_mrl, color = "BRMS"), alpha = 0.6) +
geom_abline(intercept = 0, slope = 1, color = "black", linetype = "dashed") +
scale_color_manual(values = c("AFT" = "red", "BRMS" = "blue")) +
labs(title = "MRL Prediction Comparison: AFT vs BRMS",
x = "True MRL", y = "Predicted MRL",
subtitle = paste("AFT RMSE:", round(res$aft$rmses$mrl, 3),
"| BRMS RMSE:", round(res$brms$rmses$mrl, 3))) +
theme_minimal()
print(p2)
# Bioage
bioage_comparison = data.frame(
true_b_age = brms_result$test_data$b_age,
aft_b_age = aft_result$predicted_bioage,
brms_b_age = brms_result$predicted_bioage
)
bioage_comparison = bioage_comparison[complete.cases(bioage_comparison), ]
p3 = ggplot(bioage_comparison, aes(x = true_b_age)) +
geom_point(aes(y = aft_b_age, color = "AFT"), alpha = 0.6) +
geom_point(aes(y = brms_b_age, color = "BRMS"), alpha = 0.6) +
geom_abline(intercept = 0, slope = 1, color = "black", linetype = "dashed") +
scale_color_manual(values = c("AFT" = "red", "BRMS" = "blue")) +
labs(title = "Biological Age Prediction Comparison: AFT vs BRMS",
x = "True Biological Age", y = "Predicted Biological Age",
subtitle = paste("AFT RMSE:", round(aft_result$rmses$bioage, 3),
"| BRMS RMSE:", round(brms_result$rmses$bioage, 3))) +
theme_minimal()
print(p3)
cat("\n**** SUMMARY COMPARISON ****\n")
cat("Beta RMSE - AFT:", round(res$aft$rmses$coef, 4), "\n")
cat("Beta RMSE - BRMS:", round(res$brms$rmses$coef, 4), "\n")
cat("MRL RMSE - AFT:", round(res$aft$rmses$mrl, 4), "\n")
cat("MRL RMSE - BRMS:", round(res$brms$rmses$mrl, 4), "\n")
cat("BioAge RMSE - AFT", round(res$aft$rmses$bioage, 4), "\n")
cat("BioAge RMSE - BRMS", round(res$brms$rmses$bioage, 4), "\n")
}
return(res)
}6.2 Method: Weibull; p = 20; \(n \in \{250, 500\}\); \(r_{nz} \in \{ 1\}\)
6.2.1 n = 250
source("R/generate_data.R")
debug_methods(n = 250, p = 20, g = 5, method = "weibull", seed = 40, t_SNR = 4, perf_A = T, nzg = 1) %>%
system.time()## **** SETTINGS ****
## n = 250 p = 20 g = 5 method = weibull
##
## Generating data...
## [1] "Generating βs..."
## Mean of betas pre-scaling: 0.286206588904146
## Lower range of beta values pre-scaling: -3.42901337506225
## Upper range of beta values pre-scaling: 4.47580730241499## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## ℹ Please use `as_tibble()` instead.
## ℹ The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
## `.name_repair` is omitted as of tibble 2.0.0.
## ℹ Using compatibility `.name_repair`.
## ℹ The deprecated feature was likely used in the tibble package.
## Please report the issue at <https://github.com/tidyverse/tibble/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Method: weibull
## Target SNR: 4
## Epsilon variance: 0.0165
## Mean of betas post-scaling: 2.25514051876985e-18
## Lower range of beta values post-scaling: -0.0422802473570912
## Upper range of beta values post-scaling: 0.0476788336122017
## Achieved SNR: 4
## [1] "Generating X..."
## [1] "Done!"
## [1] "Generating Weibull population lifetable..."
## Range of death: 27.7661009696673
## Range of death: 143.007815466488
## [1] "Done!"
## [1] "Generating data from Weibull lifetable..."
## Nr of valid indices: 2179
## Range of linpred values: -0.5733417 0.455935
## [1] "Done!"
## Data generated successfully!
## Dataset size: 500 x 29
##
## **** TESTING AFT ****
## [1] "----------"
## shape_est: 11.90365
## scale_est: 78.61759
## [1] "----------"
## AFT SUCCESS!
## AFT RMSEs:
## $coef
## [1] 0.01134541
##
## $mrl
## [1] 3.751361
##
## $bioage
## [1] 6.265548
##
## **** TESTING BRMS ****## Start sampling## Warning in .fun(data = .x1, seed = .x2, init = .x3, adapt_delta = .x4,
## max_treedepth = .x5, : 'stepsize' is deprecated. Please use 'step_size'
## instead.## Running MCMC with 2 parallel chains, with 2 thread(s) per chain...
##
## Chain 1 Iteration: 1 / 6000 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 6000 [ 0%] (Warmup)
## Chain 2 Iteration: 100 / 6000 [ 1%] (Warmup)
## Chain 1 Iteration: 100 / 6000 [ 1%] (Warmup)
## Chain 2 Iteration: 200 / 6000 [ 3%] (Warmup)
## Chain 2 Iteration: 300 / 6000 [ 5%] (Warmup)
## Chain 1 Iteration: 200 / 6000 [ 3%] (Warmup)
## Chain 1 Iteration: 300 / 6000 [ 5%] (Warmup)
## Chain 2 Iteration: 400 / 6000 [ 6%] (Warmup)
## Chain 1 Iteration: 400 / 6000 [ 6%] (Warmup)
## Chain 2 Iteration: 500 / 6000 [ 8%] (Warmup)
## Chain 1 Iteration: 500 / 6000 [ 8%] (Warmup)
## Chain 2 Iteration: 600 / 6000 [ 10%] (Warmup)
## Chain 1 Iteration: 600 / 6000 [ 10%] (Warmup)
## Chain 2 Iteration: 700 / 6000 [ 11%] (Warmup)
## Chain 1 Iteration: 700 / 6000 [ 11%] (Warmup)
## Chain 1 Iteration: 800 / 6000 [ 13%] (Warmup)
## Chain 2 Iteration: 800 / 6000 [ 13%] (Warmup)
## Chain 1 Iteration: 900 / 6000 [ 15%] (Warmup)
## Chain 2 Iteration: 900 / 6000 [ 15%] (Warmup)
## Chain 1 Iteration: 1000 / 6000 [ 16%] (Warmup)
## Chain 2 Iteration: 1000 / 6000 [ 16%] (Warmup)
## Chain 1 Iteration: 1100 / 6000 [ 18%] (Warmup)
## Chain 2 Iteration: 1100 / 6000 [ 18%] (Warmup)
## Chain 1 Iteration: 1200 / 6000 [ 20%] (Warmup)
## Chain 2 Iteration: 1200 / 6000 [ 20%] (Warmup)
## Chain 1 Iteration: 1300 / 6000 [ 21%] (Warmup)
## Chain 2 Iteration: 1300 / 6000 [ 21%] (Warmup)
## Chain 1 Iteration: 1400 / 6000 [ 23%] (Warmup)
## Chain 2 Iteration: 1400 / 6000 [ 23%] (Warmup)
## Chain 1 Iteration: 1500 / 6000 [ 25%] (Warmup)
## Chain 2 Iteration: 1500 / 6000 [ 25%] (Warmup)
## Chain 1 Iteration: 1600 / 6000 [ 26%] (Warmup)
## Chain 2 Iteration: 1600 / 6000 [ 26%] (Warmup)
## Chain 1 Iteration: 1700 / 6000 [ 28%] (Warmup)
## Chain 2 Iteration: 1700 / 6000 [ 28%] (Warmup)
## Chain 1 Iteration: 1800 / 6000 [ 30%] (Warmup)
## Chain 2 Iteration: 1800 / 6000 [ 30%] (Warmup)
## Chain 1 Iteration: 1900 / 6000 [ 31%] (Warmup)
## Chain 2 Iteration: 1900 / 6000 [ 31%] (Warmup)
## Chain 1 Iteration: 2000 / 6000 [ 33%] (Warmup)
## Chain 2 Iteration: 2000 / 6000 [ 33%] (Warmup)
## Chain 1 Iteration: 2100 / 6000 [ 35%] (Warmup)
## Chain 2 Iteration: 2100 / 6000 [ 35%] (Warmup)
## Chain 1 Iteration: 2200 / 6000 [ 36%] (Warmup)
## Chain 2 Iteration: 2200 / 6000 [ 36%] (Warmup)
## Chain 1 Iteration: 2300 / 6000 [ 38%] (Warmup)
## Chain 2 Iteration: 2300 / 6000 [ 38%] (Warmup)
## Chain 1 Iteration: 2400 / 6000 [ 40%] (Warmup)
## Chain 2 Iteration: 2400 / 6000 [ 40%] (Warmup)
## Chain 1 Iteration: 2500 / 6000 [ 41%] (Warmup)
## Chain 2 Iteration: 2500 / 6000 [ 41%] (Warmup)
## Chain 1 Iteration: 2600 / 6000 [ 43%] (Warmup)
## Chain 2 Iteration: 2600 / 6000 [ 43%] (Warmup)
## Chain 1 Iteration: 2700 / 6000 [ 45%] (Warmup)
## Chain 2 Iteration: 2700 / 6000 [ 45%] (Warmup)
## Chain 1 Iteration: 2800 / 6000 [ 46%] (Warmup)
## Chain 2 Iteration: 2800 / 6000 [ 46%] (Warmup)
## Chain 2 Iteration: 2900 / 6000 [ 48%] (Warmup)
## Chain 1 Iteration: 2900 / 6000 [ 48%] (Warmup)
## Chain 1 Iteration: 3000 / 6000 [ 50%] (Warmup)
## Chain 1 Iteration: 3001 / 6000 [ 50%] (Sampling)
## Chain 2 Iteration: 3000 / 6000 [ 50%] (Warmup)
## Chain 2 Iteration: 3001 / 6000 [ 50%] (Sampling)
## Chain 1 Iteration: 3100 / 6000 [ 51%] (Sampling)
## Chain 2 Iteration: 3100 / 6000 [ 51%] (Sampling)
## Chain 1 Iteration: 3200 / 6000 [ 53%] (Sampling)
## Chain 2 Iteration: 3200 / 6000 [ 53%] (Sampling)
## Chain 1 Iteration: 3300 / 6000 [ 55%] (Sampling)
## Chain 2 Iteration: 3300 / 6000 [ 55%] (Sampling)
## Chain 1 Iteration: 3400 / 6000 [ 56%] (Sampling)
## Chain 2 Iteration: 3400 / 6000 [ 56%] (Sampling)
## Chain 1 Iteration: 3500 / 6000 [ 58%] (Sampling)
## Chain 2 Iteration: 3500 / 6000 [ 58%] (Sampling)
## Chain 1 Iteration: 3600 / 6000 [ 60%] (Sampling)
## Chain 2 Iteration: 3600 / 6000 [ 60%] (Sampling)
## Chain 1 Iteration: 3700 / 6000 [ 61%] (Sampling)
## Chain 2 Iteration: 3700 / 6000 [ 61%] (Sampling)
## Chain 1 Iteration: 3800 / 6000 [ 63%] (Sampling)
## Chain 2 Iteration: 3800 / 6000 [ 63%] (Sampling)
## Chain 1 Iteration: 3900 / 6000 [ 65%] (Sampling)
## Chain 2 Iteration: 3900 / 6000 [ 65%] (Sampling)
## Chain 1 Iteration: 4000 / 6000 [ 66%] (Sampling)
## Chain 2 Iteration: 4000 / 6000 [ 66%] (Sampling)
## Chain 1 Iteration: 4100 / 6000 [ 68%] (Sampling)
## Chain 2 Iteration: 4100 / 6000 [ 68%] (Sampling)
## Chain 1 Iteration: 4200 / 6000 [ 70%] (Sampling)
## Chain 2 Iteration: 4200 / 6000 [ 70%] (Sampling)
## Chain 1 Iteration: 4300 / 6000 [ 71%] (Sampling)
## Chain 2 Iteration: 4300 / 6000 [ 71%] (Sampling)
## Chain 1 Iteration: 4400 / 6000 [ 73%] (Sampling)
## Chain 2 Iteration: 4400 / 6000 [ 73%] (Sampling)
## Chain 1 Iteration: 4500 / 6000 [ 75%] (Sampling)
## Chain 2 Iteration: 4500 / 6000 [ 75%] (Sampling)
## Chain 1 Iteration: 4600 / 6000 [ 76%] (Sampling)
## Chain 2 Iteration: 4600 / 6000 [ 76%] (Sampling)
## Chain 1 Iteration: 4700 / 6000 [ 78%] (Sampling)
## Chain 2 Iteration: 4700 / 6000 [ 78%] (Sampling)
## Chain 1 Iteration: 4800 / 6000 [ 80%] (Sampling)
## Chain 2 Iteration: 4800 / 6000 [ 80%] (Sampling)
## Chain 1 Iteration: 4900 / 6000 [ 81%] (Sampling)
## Chain 2 Iteration: 4900 / 6000 [ 81%] (Sampling)
## Chain 1 Iteration: 5000 / 6000 [ 83%] (Sampling)
## Chain 2 Iteration: 5000 / 6000 [ 83%] (Sampling)
## Chain 1 Iteration: 5100 / 6000 [ 85%] (Sampling)
## Chain 2 Iteration: 5100 / 6000 [ 85%] (Sampling)
## Chain 1 Iteration: 5200 / 6000 [ 86%] (Sampling)
## Chain 2 Iteration: 5200 / 6000 [ 86%] (Sampling)
## Chain 1 Iteration: 5300 / 6000 [ 88%] (Sampling)
## Chain 2 Iteration: 5300 / 6000 [ 88%] (Sampling)
## Chain 1 Iteration: 5400 / 6000 [ 90%] (Sampling)
## Chain 2 Iteration: 5400 / 6000 [ 90%] (Sampling)
## Chain 1 Iteration: 5500 / 6000 [ 91%] (Sampling)
## Chain 2 Iteration: 5500 / 6000 [ 91%] (Sampling)
## Chain 1 Iteration: 5600 / 6000 [ 93%] (Sampling)
## Chain 2 Iteration: 5600 / 6000 [ 93%] (Sampling)
## Chain 1 Iteration: 5700 / 6000 [ 95%] (Sampling)
## Chain 2 Iteration: 5700 / 6000 [ 95%] (Sampling)
## Chain 1 Iteration: 5800 / 6000 [ 96%] (Sampling)
## Chain 2 Iteration: 5800 / 6000 [ 96%] (Sampling)
## Chain 1 Iteration: 5900 / 6000 [ 98%] (Sampling)
## Chain 2 Iteration: 5900 / 6000 [ 98%] (Sampling)
## Chain 1 Iteration: 6000 / 6000 [100%] (Sampling)
## Chain 1 finished in 31.8 seconds.
## Chain 2 Iteration: 6000 / 6000 [100%] (Sampling)
## Chain 2 finished in 31.9 seconds.
##
## Both chains finished successfully.
## Mean chain execution time: 31.9 seconds.
## Total execution time: 32.0 seconds.## Warning: 19 of 6000 (0.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.## Loading required package: rstan## Warning: package 'rstan' was built under R version 4.3.3## Loading required package: StanHeaders## Warning: package 'StanHeaders' was built under R version 4.3.3##
## rstan version 2.32.7 (Stan version 2.32.2)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)##
## Attaching package: 'rstan'## The following object is masked from 'package:tidyr':
##
## extract## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.## Shape estimate: 9.9499
## Scale estimate 75.63284
## BRMS SUCCESS!
## BRMS RMSEs:
## $coef
## [1] 0.01258627
##
## $mrl
## [1] 5.790734
##
## $bioage
## [1] 8.876742
##
## **** COMPARISON ****
##
## **** SUMMARY COMPARISON ****
## Beta RMSE - AFT: 0.0113
## Beta RMSE - BRMS: 0.0126
## MRL RMSE - AFT: 3.7514
## MRL RMSE - BRMS: 5.7907
## BioAge RMSE - AFT 6.2655
## BioAge RMSE - BRMS 8.8767## user system elapsed
## 117.908 11.205 45.8296.2.2 n = 500
debug_methods(n = 500, p = 20, g = 5, method = "weibull", seed = 40, t_SNR = 4, perf_A = T, nzg = 1) %>%
system.time()## **** SETTINGS ****
## n = 500 p = 20 g = 5 method = weibull
##
## Generating data...
## [1] "Generating βs..."
## Mean of betas pre-scaling: 0.286206588904146
## Lower range of beta values pre-scaling: -3.42901337506225
## Upper range of beta values pre-scaling: 4.47580730241499
## Method: weibull
## Target SNR: 4
## Epsilon variance: 0.0165
## Mean of betas post-scaling: 2.25514051876985e-18
## Lower range of beta values post-scaling: -0.0422802473570912
## Upper range of beta values post-scaling: 0.0476788336122017
## Achieved SNR: 4
## [1] "Generating X..."
## [1] "Done!"
## [1] "Generating Weibull population lifetable..."
## Range of death: 27.7661009696673
## Range of death: 143.007815466488
## [1] "Done!"
## [1] "Generating data from Weibull lifetable..."
## Nr of valid indices: 4407
## Range of linpred values: -0.5280914 0.5125125
## [1] "Done!"
## Data generated successfully!
## Dataset size: 1000 x 29
##
## **** TESTING AFT ****
## [1] "----------"
## shape_est: 10.06515
## scale_est: 80.82797
## [1] "----------"
## AFT SUCCESS!
## AFT RMSEs:
## $coef
## [1] 0.009912031
##
## $mrl
## [1] 3.294007
##
## $bioage
## [1] 6.440819
##
## **** TESTING BRMS ****## Start sampling## Warning in .fun(data = .x1, seed = .x2, init = .x3, adapt_delta = .x4,
## max_treedepth = .x5, : 'stepsize' is deprecated. Please use 'step_size'
## instead.## Running MCMC with 2 parallel chains, with 2 thread(s) per chain...
##
## Chain 1 Iteration: 1 / 6000 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 6000 [ 0%] (Warmup)
## Chain 1 Iteration: 100 / 6000 [ 1%] (Warmup)
## Chain 2 Iteration: 100 / 6000 [ 1%] (Warmup)
## Chain 1 Iteration: 200 / 6000 [ 3%] (Warmup)
## Chain 2 Iteration: 200 / 6000 [ 3%] (Warmup)
## Chain 1 Iteration: 300 / 6000 [ 5%] (Warmup)
## Chain 2 Iteration: 300 / 6000 [ 5%] (Warmup)
## Chain 1 Iteration: 400 / 6000 [ 6%] (Warmup)
## Chain 2 Iteration: 400 / 6000 [ 6%] (Warmup)
## Chain 2 Iteration: 500 / 6000 [ 8%] (Warmup)
## Chain 1 Iteration: 500 / 6000 [ 8%] (Warmup)
## Chain 1 Iteration: 600 / 6000 [ 10%] (Warmup)
## Chain 2 Iteration: 600 / 6000 [ 10%] (Warmup)
## Chain 2 Iteration: 700 / 6000 [ 11%] (Warmup)
## Chain 1 Iteration: 700 / 6000 [ 11%] (Warmup)
## Chain 2 Iteration: 800 / 6000 [ 13%] (Warmup)
## Chain 1 Iteration: 800 / 6000 [ 13%] (Warmup)
## Chain 2 Iteration: 900 / 6000 [ 15%] (Warmup)
## Chain 1 Iteration: 900 / 6000 [ 15%] (Warmup)
## Chain 2 Iteration: 1000 / 6000 [ 16%] (Warmup)
## Chain 1 Iteration: 1000 / 6000 [ 16%] (Warmup)
## Chain 2 Iteration: 1100 / 6000 [ 18%] (Warmup)
## Chain 1 Iteration: 1100 / 6000 [ 18%] (Warmup)
## Chain 2 Iteration: 1200 / 6000 [ 20%] (Warmup)
## Chain 1 Iteration: 1200 / 6000 [ 20%] (Warmup)
## Chain 2 Iteration: 1300 / 6000 [ 21%] (Warmup)
## Chain 1 Iteration: 1300 / 6000 [ 21%] (Warmup)
## Chain 2 Iteration: 1400 / 6000 [ 23%] (Warmup)
## Chain 1 Iteration: 1400 / 6000 [ 23%] (Warmup)
## Chain 2 Iteration: 1500 / 6000 [ 25%] (Warmup)
## Chain 1 Iteration: 1500 / 6000 [ 25%] (Warmup)
## Chain 2 Iteration: 1600 / 6000 [ 26%] (Warmup)
## Chain 1 Iteration: 1600 / 6000 [ 26%] (Warmup)
## Chain 2 Iteration: 1700 / 6000 [ 28%] (Warmup)
## Chain 1 Iteration: 1700 / 6000 [ 28%] (Warmup)
## Chain 2 Iteration: 1800 / 6000 [ 30%] (Warmup)
## Chain 1 Iteration: 1800 / 6000 [ 30%] (Warmup)
## Chain 2 Iteration: 1900 / 6000 [ 31%] (Warmup)
## Chain 1 Iteration: 1900 / 6000 [ 31%] (Warmup)
## Chain 2 Iteration: 2000 / 6000 [ 33%] (Warmup)
## Chain 1 Iteration: 2000 / 6000 [ 33%] (Warmup)
## Chain 2 Iteration: 2100 / 6000 [ 35%] (Warmup)
## Chain 1 Iteration: 2100 / 6000 [ 35%] (Warmup)
## Chain 2 Iteration: 2200 / 6000 [ 36%] (Warmup)
## Chain 1 Iteration: 2200 / 6000 [ 36%] (Warmup)
## Chain 2 Iteration: 2300 / 6000 [ 38%] (Warmup)
## Chain 1 Iteration: 2300 / 6000 [ 38%] (Warmup)
## Chain 2 Iteration: 2400 / 6000 [ 40%] (Warmup)
## Chain 1 Iteration: 2400 / 6000 [ 40%] (Warmup)
## Chain 2 Iteration: 2500 / 6000 [ 41%] (Warmup)
## Chain 1 Iteration: 2500 / 6000 [ 41%] (Warmup)
## Chain 2 Iteration: 2600 / 6000 [ 43%] (Warmup)
## Chain 1 Iteration: 2600 / 6000 [ 43%] (Warmup)
## Chain 2 Iteration: 2700 / 6000 [ 45%] (Warmup)
## Chain 1 Iteration: 2700 / 6000 [ 45%] (Warmup)
## Chain 2 Iteration: 2800 / 6000 [ 46%] (Warmup)
## Chain 1 Iteration: 2800 / 6000 [ 46%] (Warmup)
## Chain 2 Iteration: 2900 / 6000 [ 48%] (Warmup)
## Chain 1 Iteration: 2900 / 6000 [ 48%] (Warmup)
## Chain 2 Iteration: 3000 / 6000 [ 50%] (Warmup)
## Chain 2 Iteration: 3001 / 6000 [ 50%] (Sampling)
## Chain 1 Iteration: 3000 / 6000 [ 50%] (Warmup)
## Chain 1 Iteration: 3001 / 6000 [ 50%] (Sampling)
## Chain 2 Iteration: 3100 / 6000 [ 51%] (Sampling)
## Chain 1 Iteration: 3100 / 6000 [ 51%] (Sampling)
## Chain 2 Iteration: 3200 / 6000 [ 53%] (Sampling)
## Chain 1 Iteration: 3200 / 6000 [ 53%] (Sampling)
## Chain 2 Iteration: 3300 / 6000 [ 55%] (Sampling)
## Chain 1 Iteration: 3300 / 6000 [ 55%] (Sampling)
## Chain 2 Iteration: 3400 / 6000 [ 56%] (Sampling)
## Chain 1 Iteration: 3400 / 6000 [ 56%] (Sampling)
## Chain 2 Iteration: 3500 / 6000 [ 58%] (Sampling)
## Chain 1 Iteration: 3500 / 6000 [ 58%] (Sampling)
## Chain 2 Iteration: 3600 / 6000 [ 60%] (Sampling)
## Chain 1 Iteration: 3600 / 6000 [ 60%] (Sampling)
## Chain 2 Iteration: 3700 / 6000 [ 61%] (Sampling)
## Chain 1 Iteration: 3700 / 6000 [ 61%] (Sampling)
## Chain 2 Iteration: 3800 / 6000 [ 63%] (Sampling)
## Chain 1 Iteration: 3800 / 6000 [ 63%] (Sampling)
## Chain 2 Iteration: 3900 / 6000 [ 65%] (Sampling)
## Chain 1 Iteration: 3900 / 6000 [ 65%] (Sampling)
## Chain 2 Iteration: 4000 / 6000 [ 66%] (Sampling)
## Chain 1 Iteration: 4000 / 6000 [ 66%] (Sampling)
## Chain 2 Iteration: 4100 / 6000 [ 68%] (Sampling)
## Chain 2 Iteration: 4200 / 6000 [ 70%] (Sampling)
## Chain 1 Iteration: 4100 / 6000 [ 68%] (Sampling)
## Chain 2 Iteration: 4300 / 6000 [ 71%] (Sampling)
## Chain 1 Iteration: 4200 / 6000 [ 70%] (Sampling)
## Chain 2 Iteration: 4400 / 6000 [ 73%] (Sampling)
## Chain 1 Iteration: 4300 / 6000 [ 71%] (Sampling)
## Chain 2 Iteration: 4500 / 6000 [ 75%] (Sampling)
## Chain 1 Iteration: 4400 / 6000 [ 73%] (Sampling)
## Chain 2 Iteration: 4600 / 6000 [ 76%] (Sampling)
## Chain 1 Iteration: 4500 / 6000 [ 75%] (Sampling)
## Chain 2 Iteration: 4700 / 6000 [ 78%] (Sampling)
## Chain 1 Iteration: 4600 / 6000 [ 76%] (Sampling)
## Chain 2 Iteration: 4800 / 6000 [ 80%] (Sampling)
## Chain 1 Iteration: 4700 / 6000 [ 78%] (Sampling)
## Chain 2 Iteration: 4900 / 6000 [ 81%] (Sampling)
## Chain 1 Iteration: 4800 / 6000 [ 80%] (Sampling)
## Chain 2 Iteration: 5000 / 6000 [ 83%] (Sampling)
## Chain 1 Iteration: 4900 / 6000 [ 81%] (Sampling)
## Chain 2 Iteration: 5100 / 6000 [ 85%] (Sampling)
## Chain 1 Iteration: 5000 / 6000 [ 83%] (Sampling)
## Chain 2 Iteration: 5200 / 6000 [ 86%] (Sampling)
## Chain 1 Iteration: 5100 / 6000 [ 85%] (Sampling)
## Chain 2 Iteration: 5300 / 6000 [ 88%] (Sampling)
## Chain 1 Iteration: 5200 / 6000 [ 86%] (Sampling)
## Chain 2 Iteration: 5400 / 6000 [ 90%] (Sampling)
## Chain 1 Iteration: 5300 / 6000 [ 88%] (Sampling)
## Chain 2 Iteration: 5500 / 6000 [ 91%] (Sampling)
## Chain 1 Iteration: 5400 / 6000 [ 90%] (Sampling)
## Chain 2 Iteration: 5600 / 6000 [ 93%] (Sampling)
## Chain 1 Iteration: 5500 / 6000 [ 91%] (Sampling)
## Chain 2 Iteration: 5700 / 6000 [ 95%] (Sampling)
## Chain 1 Iteration: 5600 / 6000 [ 93%] (Sampling)
## Chain 2 Iteration: 5800 / 6000 [ 96%] (Sampling)
## Chain 1 Iteration: 5700 / 6000 [ 95%] (Sampling)
## Chain 2 Iteration: 5900 / 6000 [ 98%] (Sampling)
## Chain 1 Iteration: 5800 / 6000 [ 96%] (Sampling)
## Chain 2 Iteration: 6000 / 6000 [100%] (Sampling)
## Chain 2 finished in 79.4 seconds.
## Chain 1 Iteration: 5900 / 6000 [ 98%] (Sampling)
## Chain 1 Iteration: 6000 / 6000 [100%] (Sampling)
## Chain 1 finished in 80.8 seconds.
##
## Both chains finished successfully.
## Mean chain execution time: 80.1 seconds.
## Total execution time: 80.9 seconds.## Warning: 42 of 6000 (1.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.## Shape estimate: 9.378186
## Scale estimate 76.9985
## BRMS SUCCESS!
## BRMS RMSEs:
## $coef
## [1] 0.01035903
##
## $mrl
## [1] 4.586242
##
## $bioage
## [1] 8.699857
##
## **** COMPARISON ****
##
## **** SUMMARY COMPARISON ****
## Beta RMSE - AFT: 0.0099
## Beta RMSE - BRMS: 0.0104
## MRL RMSE - AFT: 3.294
## MRL RMSE - BRMS: 4.5862
## BioAge RMSE - AFT 6.4408
## BioAge RMSE - BRMS 8.6999## user system elapsed
## 264.076 22.496 85.0826.3 Method: Weibull; p = 100; \(n \in \{250, 500\}\); \(r_{nz} \in \{1\}\)
6.3.1 n = 250
debug_methods(n = 250, p = 100, g = 5, method = "weibull", seed = 40, t_SNR = 4, perf_A = T, nzg = 1) %>%
system.time()## **** SETTINGS ****
## n = 250 p = 100 g = 5 method = weibull
##
## Generating data...
## [1] "Generating βs..."
## Mean of betas pre-scaling: 0.249922847059995
## Lower range of beta values pre-scaling: -3.52945436691976
## Upper range of beta values pre-scaling: 4.83157475138873
## Method: weibull
## Target SNR: 4
## Epsilon variance: 0.0165
## Mean of betas post-scaling: 2.42861286636753e-19
## Lower range of beta values post-scaling: -0.0112344086238269
## Upper range of beta values post-scaling: 0.0136192146883276
## Achieved SNR: 4
## [1] "Generating X..."
## [1] "Done!"
## [1] "Generating Weibull population lifetable..."
## Range of death: 29.1486009954669
## Range of death: 113.639483758283
## [1] "Done!"
## [1] "Generating data from Weibull lifetable..."
## Nr of valid indices: 2240
## Range of linpred values: -0.2851693 0.2659974
## [1] "Done!"
## Data generated successfully!
## Dataset size: 500 x 109
##
## **** TESTING AFT ****## AFT model error: Singular hessian; suspicious variable No. TRUE:
## =
## Try running with fixed shape## [1] "Fit failed"
## AFT FAILED!
##
## **** TESTING BRMS ****## Start sampling## Warning in .fun(data = .x1, seed = .x2, init = .x3, adapt_delta = .x4,
## max_treedepth = .x5, : 'stepsize' is deprecated. Please use 'step_size'
## instead.## Running MCMC with 2 parallel chains, with 2 thread(s) per chain...
##
## Chain 1 Iteration: 1 / 6000 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 6000 [ 0%] (Warmup)
## Chain 1 Iteration: 100 / 6000 [ 1%] (Warmup)
## Chain 2 Iteration: 100 / 6000 [ 1%] (Warmup)
## Chain 2 Iteration: 200 / 6000 [ 3%] (Warmup)
## Chain 1 Iteration: 200 / 6000 [ 3%] (Warmup)
## Chain 2 Iteration: 300 / 6000 [ 5%] (Warmup)
## Chain 1 Iteration: 300 / 6000 [ 5%] (Warmup)
## Chain 2 Iteration: 400 / 6000 [ 6%] (Warmup)
## Chain 1 Iteration: 400 / 6000 [ 6%] (Warmup)
## Chain 2 Iteration: 500 / 6000 [ 8%] (Warmup)
## Chain 1 Iteration: 500 / 6000 [ 8%] (Warmup)
## Chain 2 Iteration: 600 / 6000 [ 10%] (Warmup)
## Chain 1 Iteration: 600 / 6000 [ 10%] (Warmup)
## Chain 2 Iteration: 700 / 6000 [ 11%] (Warmup)
## Chain 2 Iteration: 800 / 6000 [ 13%] (Warmup)
## Chain 1 Iteration: 700 / 6000 [ 11%] (Warmup)
## Chain 1 Iteration: 800 / 6000 [ 13%] (Warmup)
## Chain 2 Iteration: 900 / 6000 [ 15%] (Warmup)
## Chain 2 Iteration: 1000 / 6000 [ 16%] (Warmup)
## Chain 1 Iteration: 900 / 6000 [ 15%] (Warmup)
## Chain 2 Iteration: 1100 / 6000 [ 18%] (Warmup)
## Chain 1 Iteration: 1000 / 6000 [ 16%] (Warmup)
## Chain 2 Iteration: 1200 / 6000 [ 20%] (Warmup)
## Chain 1 Iteration: 1100 / 6000 [ 18%] (Warmup)
## Chain 2 Iteration: 1300 / 6000 [ 21%] (Warmup)
## Chain 2 Iteration: 1400 / 6000 [ 23%] (Warmup)
## Chain 1 Iteration: 1200 / 6000 [ 20%] (Warmup)
## Chain 1 Iteration: 1300 / 6000 [ 21%] (Warmup)
## Chain 2 Iteration: 1500 / 6000 [ 25%] (Warmup)
## Chain 1 Iteration: 1400 / 6000 [ 23%] (Warmup)
## Chain 2 Iteration: 1600 / 6000 [ 26%] (Warmup)
## Chain 1 Iteration: 1500 / 6000 [ 25%] (Warmup)
## Chain 2 Iteration: 1700 / 6000 [ 28%] (Warmup)
## Chain 1 Iteration: 1600 / 6000 [ 26%] (Warmup)
## Chain 2 Iteration: 1800 / 6000 [ 30%] (Warmup)
## Chain 1 Iteration: 1700 / 6000 [ 28%] (Warmup)
## Chain 2 Iteration: 1900 / 6000 [ 31%] (Warmup)
## Chain 1 Iteration: 1800 / 6000 [ 30%] (Warmup)
## Chain 2 Iteration: 2000 / 6000 [ 33%] (Warmup)
## Chain 1 Iteration: 1900 / 6000 [ 31%] (Warmup)
## Chain 1 Iteration: 2000 / 6000 [ 33%] (Warmup)
## Chain 2 Iteration: 2100 / 6000 [ 35%] (Warmup)
## Chain 1 Iteration: 2100 / 6000 [ 35%] (Warmup)
## Chain 2 Iteration: 2200 / 6000 [ 36%] (Warmup)
## Chain 1 Iteration: 2200 / 6000 [ 36%] (Warmup)
## Chain 2 Iteration: 2300 / 6000 [ 38%] (Warmup)
## Chain 2 Iteration: 2400 / 6000 [ 40%] (Warmup)
## Chain 1 Iteration: 2300 / 6000 [ 38%] (Warmup)
## Chain 2 Iteration: 2500 / 6000 [ 41%] (Warmup)
## Chain 1 Iteration: 2400 / 6000 [ 40%] (Warmup)
## Chain 2 Iteration: 2600 / 6000 [ 43%] (Warmup)
## Chain 1 Iteration: 2500 / 6000 [ 41%] (Warmup)
## Chain 2 Iteration: 2700 / 6000 [ 45%] (Warmup)
## Chain 1 Iteration: 2600 / 6000 [ 43%] (Warmup)
## Chain 2 Iteration: 2800 / 6000 [ 46%] (Warmup)
## Chain 1 Iteration: 2700 / 6000 [ 45%] (Warmup)
## Chain 2 Iteration: 2900 / 6000 [ 48%] (Warmup)
## Chain 1 Iteration: 2800 / 6000 [ 46%] (Warmup)
## Chain 2 Iteration: 3000 / 6000 [ 50%] (Warmup)
## Chain 2 Iteration: 3001 / 6000 [ 50%] (Sampling)
## Chain 1 Iteration: 2900 / 6000 [ 48%] (Warmup)
## Chain 1 Iteration: 3000 / 6000 [ 50%] (Warmup)
## Chain 2 Iteration: 3100 / 6000 [ 51%] (Sampling)
## Chain 1 Iteration: 3001 / 6000 [ 50%] (Sampling)
## Chain 1 Iteration: 3100 / 6000 [ 51%] (Sampling)
## Chain 2 Iteration: 3200 / 6000 [ 53%] (Sampling)
## Chain 1 Iteration: 3200 / 6000 [ 53%] (Sampling)
## Chain 1 Iteration: 3300 / 6000 [ 55%] (Sampling)
## Chain 2 Iteration: 3300 / 6000 [ 55%] (Sampling)
## Chain 1 Iteration: 3400 / 6000 [ 56%] (Sampling)
## Chain 1 Iteration: 3500 / 6000 [ 58%] (Sampling)
## Chain 2 Iteration: 3400 / 6000 [ 56%] (Sampling)
## Chain 1 Iteration: 3600 / 6000 [ 60%] (Sampling)
## Chain 1 Iteration: 3700 / 6000 [ 61%] (Sampling)
## Chain 2 Iteration: 3500 / 6000 [ 58%] (Sampling)
## Chain 1 Iteration: 3800 / 6000 [ 63%] (Sampling)
## Chain 1 Iteration: 3900 / 6000 [ 65%] (Sampling)
## Chain 2 Iteration: 3600 / 6000 [ 60%] (Sampling)
## Chain 1 Iteration: 4000 / 6000 [ 66%] (Sampling)
## Chain 2 Iteration: 3700 / 6000 [ 61%] (Sampling)
## Chain 1 Iteration: 4100 / 6000 [ 68%] (Sampling)
## Chain 1 Iteration: 4200 / 6000 [ 70%] (Sampling)
## Chain 2 Iteration: 3800 / 6000 [ 63%] (Sampling)
## Chain 1 Iteration: 4300 / 6000 [ 71%] (Sampling)
## Chain 1 Iteration: 4400 / 6000 [ 73%] (Sampling)
## Chain 2 Iteration: 3900 / 6000 [ 65%] (Sampling)
## Chain 1 Iteration: 4500 / 6000 [ 75%] (Sampling)
## Chain 1 Iteration: 4600 / 6000 [ 76%] (Sampling)
## Chain 2 Iteration: 4000 / 6000 [ 66%] (Sampling)
## Chain 1 Iteration: 4700 / 6000 [ 78%] (Sampling)
## Chain 1 Iteration: 4800 / 6000 [ 80%] (Sampling)
## Chain 2 Iteration: 4100 / 6000 [ 68%] (Sampling)
## Chain 1 Iteration: 4900 / 6000 [ 81%] (Sampling)
## Chain 2 Iteration: 4200 / 6000 [ 70%] (Sampling)
## Chain 1 Iteration: 5000 / 6000 [ 83%] (Sampling)
## Chain 1 Iteration: 5100 / 6000 [ 85%] (Sampling)
## Chain 2 Iteration: 4300 / 6000 [ 71%] (Sampling)
## Chain 1 Iteration: 5200 / 6000 [ 86%] (Sampling)
## Chain 1 Iteration: 5300 / 6000 [ 88%] (Sampling)
## Chain 2 Iteration: 4400 / 6000 [ 73%] (Sampling)
## Chain 1 Iteration: 5400 / 6000 [ 90%] (Sampling)
## Chain 1 Iteration: 5500 / 6000 [ 91%] (Sampling)
## Chain 2 Iteration: 4500 / 6000 [ 75%] (Sampling)
## Chain 1 Iteration: 5600 / 6000 [ 93%] (Sampling)
## Chain 1 Iteration: 5700 / 6000 [ 95%] (Sampling)
## Chain 2 Iteration: 4600 / 6000 [ 76%] (Sampling)
## Chain 1 Iteration: 5800 / 6000 [ 96%] (Sampling)
## Chain 1 Iteration: 5900 / 6000 [ 98%] (Sampling)
## Chain 2 Iteration: 4700 / 6000 [ 78%] (Sampling)
## Chain 1 Iteration: 6000 / 6000 [100%] (Sampling)
## Chain 1 finished in 32.8 seconds.
## Chain 2 Iteration: 4800 / 6000 [ 80%] (Sampling)
## Chain 2 Iteration: 4900 / 6000 [ 81%] (Sampling)
## Chain 2 Iteration: 5000 / 6000 [ 83%] (Sampling)
## Chain 2 Iteration: 5100 / 6000 [ 85%] (Sampling)
## Chain 2 Iteration: 5200 / 6000 [ 86%] (Sampling)
## Chain 2 Iteration: 5300 / 6000 [ 88%] (Sampling)
## Chain 2 Iteration: 5400 / 6000 [ 90%] (Sampling)
## Chain 2 Iteration: 5500 / 6000 [ 91%] (Sampling)
## Chain 2 Iteration: 5600 / 6000 [ 93%] (Sampling)
## Chain 2 Iteration: 5700 / 6000 [ 95%] (Sampling)
## Chain 2 Iteration: 5800 / 6000 [ 96%] (Sampling)
## Chain 2 Iteration: 5900 / 6000 [ 98%] (Sampling)
## Chain 2 Iteration: 6000 / 6000 [100%] (Sampling)
## Chain 2 finished in 40.1 seconds.
##
## Both chains finished successfully.
## Mean chain execution time: 36.5 seconds.
## Total execution time: 40.2 seconds.## Warning: 1 of 6000 (0.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.## Shape estimate: 8.818291
## Scale estimate 76.27433
## BRMS SUCCESS!
## BRMS RMSEs:
## $coef
## [1] 0.007108537
##
## $mrl
## [1] 6.190739
##
## $bioage
## [1] 7.873618
##
## **** COMPARISON ****
##
## **** SUMMARY COMPARISON ****
## Beta RMSE - AFT: 0
## Beta RMSE - BRMS: 0.0071
## MRL RMSE - AFT: 0
## MRL RMSE - BRMS: 6.1907
## BioAge RMSE - AFT 0
## BioAge RMSE - BRMS 7.8736## user system elapsed
## 128.922 13.561 54.6996.3.2 n = 500
debug_methods(n = 500, p = 100, g = 5, method = "weibull", seed = 40, t_SNR = 4, perf_A = T, nzg = 1) %>%
system.time()## **** SETTINGS ****
## n = 500 p = 100 g = 5 method = weibull
##
## Generating data...
## [1] "Generating βs..."
## Mean of betas pre-scaling: 0.249922847059995
## Lower range of beta values pre-scaling: -3.52945436691976
## Upper range of beta values pre-scaling: 4.83157475138873
## Method: weibull
## Target SNR: 4
## Epsilon variance: 0.0165
## Mean of betas post-scaling: 2.42861286636753e-19
## Lower range of beta values post-scaling: -0.0112344086238269
## Upper range of beta values post-scaling: 0.0136192146883276
## Achieved SNR: 4
## [1] "Generating X..."
## [1] "Done!"
## [1] "Generating Weibull population lifetable..."
## Range of death: 29.1486009954669
## Range of death: 113.639483758283
## [1] "Done!"
## [1] "Generating data from Weibull lifetable..."
## Nr of valid indices: 4489
## Range of linpred values: -0.2702238 0.2766213
## [1] "Done!"
## Data generated successfully!
## Dataset size: 1000 x 109
##
## **** TESTING AFT ****## AFT model error: Singular hessian; suspicious variable No. TRUE:
## =
## Try running with fixed shape## [1] "Fit failed"
## AFT FAILED!
##
## **** TESTING BRMS ****## Start sampling## Warning in .fun(data = .x1, seed = .x2, init = .x3, adapt_delta = .x4,
## max_treedepth = .x5, : 'stepsize' is deprecated. Please use 'step_size'
## instead.## Running MCMC with 2 parallel chains, with 2 thread(s) per chain...
##
## Chain 1 Iteration: 1 / 6000 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 6000 [ 0%] (Warmup)
## Chain 1 Iteration: 100 / 6000 [ 1%] (Warmup)
## Chain 2 Iteration: 100 / 6000 [ 1%] (Warmup)
## Chain 1 Iteration: 200 / 6000 [ 3%] (Warmup)
## Chain 1 Iteration: 300 / 6000 [ 5%] (Warmup)
## Chain 2 Iteration: 200 / 6000 [ 3%] (Warmup)
## Chain 1 Iteration: 400 / 6000 [ 6%] (Warmup)
## Chain 2 Iteration: 300 / 6000 [ 5%] (Warmup)
## Chain 1 Iteration: 500 / 6000 [ 8%] (Warmup)
## Chain 2 Iteration: 400 / 6000 [ 6%] (Warmup)
## Chain 1 Iteration: 600 / 6000 [ 10%] (Warmup)
## Chain 2 Iteration: 500 / 6000 [ 8%] (Warmup)
## Chain 1 Iteration: 700 / 6000 [ 11%] (Warmup)
## Chain 2 Iteration: 600 / 6000 [ 10%] (Warmup)
## Chain 1 Iteration: 800 / 6000 [ 13%] (Warmup)
## Chain 2 Iteration: 700 / 6000 [ 11%] (Warmup)
## Chain 1 Iteration: 900 / 6000 [ 15%] (Warmup)
## Chain 2 Iteration: 800 / 6000 [ 13%] (Warmup)
## Chain 1 Iteration: 1000 / 6000 [ 16%] (Warmup)
## Chain 2 Iteration: 900 / 6000 [ 15%] (Warmup)
## Chain 1 Iteration: 1100 / 6000 [ 18%] (Warmup)
## Chain 2 Iteration: 1000 / 6000 [ 16%] (Warmup)
## Chain 1 Iteration: 1200 / 6000 [ 20%] (Warmup)
## Chain 2 Iteration: 1100 / 6000 [ 18%] (Warmup)
## Chain 1 Iteration: 1300 / 6000 [ 21%] (Warmup)
## Chain 2 Iteration: 1200 / 6000 [ 20%] (Warmup)
## Chain 1 Iteration: 1400 / 6000 [ 23%] (Warmup)
## Chain 2 Iteration: 1300 / 6000 [ 21%] (Warmup)
## Chain 1 Iteration: 1500 / 6000 [ 25%] (Warmup)
## Chain 2 Iteration: 1400 / 6000 [ 23%] (Warmup)
## Chain 1 Iteration: 1600 / 6000 [ 26%] (Warmup)
## Chain 2 Iteration: 1500 / 6000 [ 25%] (Warmup)
## Chain 1 Iteration: 1700 / 6000 [ 28%] (Warmup)
## Chain 2 Iteration: 1600 / 6000 [ 26%] (Warmup)
## Chain 1 Iteration: 1800 / 6000 [ 30%] (Warmup)
## Chain 2 Iteration: 1700 / 6000 [ 28%] (Warmup)
## Chain 1 Iteration: 1900 / 6000 [ 31%] (Warmup)
## Chain 2 Iteration: 1800 / 6000 [ 30%] (Warmup)
## Chain 1 Iteration: 2000 / 6000 [ 33%] (Warmup)
## Chain 2 Iteration: 1900 / 6000 [ 31%] (Warmup)
## Chain 1 Iteration: 2100 / 6000 [ 35%] (Warmup)
## Chain 2 Iteration: 2000 / 6000 [ 33%] (Warmup)
## Chain 1 Iteration: 2200 / 6000 [ 36%] (Warmup)
## Chain 2 Iteration: 2100 / 6000 [ 35%] (Warmup)
## Chain 1 Iteration: 2300 / 6000 [ 38%] (Warmup)
## Chain 2 Iteration: 2200 / 6000 [ 36%] (Warmup)
## Chain 1 Iteration: 2400 / 6000 [ 40%] (Warmup)
## Chain 2 Iteration: 2300 / 6000 [ 38%] (Warmup)
## Chain 1 Iteration: 2500 / 6000 [ 41%] (Warmup)
## Chain 2 Iteration: 2400 / 6000 [ 40%] (Warmup)
## Chain 1 Iteration: 2600 / 6000 [ 43%] (Warmup)
## Chain 2 Iteration: 2500 / 6000 [ 41%] (Warmup)
## Chain 1 Iteration: 2700 / 6000 [ 45%] (Warmup)
## Chain 2 Iteration: 2600 / 6000 [ 43%] (Warmup)
## Chain 1 Iteration: 2800 / 6000 [ 46%] (Warmup)
## Chain 2 Iteration: 2700 / 6000 [ 45%] (Warmup)
## Chain 1 Iteration: 2900 / 6000 [ 48%] (Warmup)
## Chain 2 Iteration: 2800 / 6000 [ 46%] (Warmup)
## Chain 1 Iteration: 3000 / 6000 [ 50%] (Warmup)
## Chain 1 Iteration: 3001 / 6000 [ 50%] (Sampling)
## Chain 2 Iteration: 2900 / 6000 [ 48%] (Warmup)
## Chain 1 Iteration: 3100 / 6000 [ 51%] (Sampling)
## Chain 2 Iteration: 3000 / 6000 [ 50%] (Warmup)
## Chain 2 Iteration: 3001 / 6000 [ 50%] (Sampling)
## Chain 1 Iteration: 3200 / 6000 [ 53%] (Sampling)
## Chain 2 Iteration: 3100 / 6000 [ 51%] (Sampling)
## Chain 1 Iteration: 3300 / 6000 [ 55%] (Sampling)
## Chain 2 Iteration: 3200 / 6000 [ 53%] (Sampling)
## Chain 1 Iteration: 3400 / 6000 [ 56%] (Sampling)
## Chain 2 Iteration: 3300 / 6000 [ 55%] (Sampling)
## Chain 1 Iteration: 3500 / 6000 [ 58%] (Sampling)
## Chain 2 Iteration: 3400 / 6000 [ 56%] (Sampling)
## Chain 1 Iteration: 3600 / 6000 [ 60%] (Sampling)
## Chain 2 Iteration: 3500 / 6000 [ 58%] (Sampling)
## Chain 1 Iteration: 3700 / 6000 [ 61%] (Sampling)
## Chain 2 Iteration: 3600 / 6000 [ 60%] (Sampling)
## Chain 1 Iteration: 3800 / 6000 [ 63%] (Sampling)
## Chain 2 Iteration: 3700 / 6000 [ 61%] (Sampling)
## Chain 1 Iteration: 3900 / 6000 [ 65%] (Sampling)
## Chain 2 Iteration: 3800 / 6000 [ 63%] (Sampling)
## Chain 1 Iteration: 4000 / 6000 [ 66%] (Sampling)
## Chain 2 Iteration: 3900 / 6000 [ 65%] (Sampling)
## Chain 1 Iteration: 4100 / 6000 [ 68%] (Sampling)
## Chain 1 Iteration: 4200 / 6000 [ 70%] (Sampling)
## Chain 2 Iteration: 4000 / 6000 [ 66%] (Sampling)
## Chain 1 Iteration: 4300 / 6000 [ 71%] (Sampling)
## Chain 2 Iteration: 4100 / 6000 [ 68%] (Sampling)
## Chain 2 Iteration: 4200 / 6000 [ 70%] (Sampling)
## Chain 1 Iteration: 4400 / 6000 [ 73%] (Sampling)
## Chain 1 Iteration: 4500 / 6000 [ 75%] (Sampling)
## Chain 2 Iteration: 4300 / 6000 [ 71%] (Sampling)
## Chain 1 Iteration: 4600 / 6000 [ 76%] (Sampling)
## Chain 2 Iteration: 4400 / 6000 [ 73%] (Sampling)
## Chain 1 Iteration: 4700 / 6000 [ 78%] (Sampling)
## Chain 2 Iteration: 4500 / 6000 [ 75%] (Sampling)
## Chain 1 Iteration: 4800 / 6000 [ 80%] (Sampling)
## Chain 2 Iteration: 4600 / 6000 [ 76%] (Sampling)
## Chain 1 Iteration: 4900 / 6000 [ 81%] (Sampling)
## Chain 2 Iteration: 4700 / 6000 [ 78%] (Sampling)
## Chain 1 Iteration: 5000 / 6000 [ 83%] (Sampling)
## Chain 2 Iteration: 4800 / 6000 [ 80%] (Sampling)
## Chain 1 Iteration: 5100 / 6000 [ 85%] (Sampling)
## Chain 2 Iteration: 4900 / 6000 [ 81%] (Sampling)
## Chain 1 Iteration: 5200 / 6000 [ 86%] (Sampling)
## Chain 2 Iteration: 5000 / 6000 [ 83%] (Sampling)
## Chain 1 Iteration: 5300 / 6000 [ 88%] (Sampling)
## Chain 2 Iteration: 5100 / 6000 [ 85%] (Sampling)
## Chain 1 Iteration: 5400 / 6000 [ 90%] (Sampling)
## Chain 2 Iteration: 5200 / 6000 [ 86%] (Sampling)
## Chain 1 Iteration: 5500 / 6000 [ 91%] (Sampling)
## Chain 2 Iteration: 5300 / 6000 [ 88%] (Sampling)
## Chain 1 Iteration: 5600 / 6000 [ 93%] (Sampling)
## Chain 2 Iteration: 5400 / 6000 [ 90%] (Sampling)
## Chain 1 Iteration: 5700 / 6000 [ 95%] (Sampling)
## Chain 2 Iteration: 5500 / 6000 [ 91%] (Sampling)
## Chain 1 Iteration: 5800 / 6000 [ 96%] (Sampling)
## Chain 2 Iteration: 5600 / 6000 [ 93%] (Sampling)
## Chain 1 Iteration: 5900 / 6000 [ 98%] (Sampling)
## Chain 2 Iteration: 5700 / 6000 [ 95%] (Sampling)
## Chain 1 Iteration: 6000 / 6000 [100%] (Sampling)
## Chain 1 finished in 92.3 seconds.
## Chain 2 Iteration: 5800 / 6000 [ 96%] (Sampling)
## Chain 2 Iteration: 5900 / 6000 [ 98%] (Sampling)
## Chain 2 Iteration: 6000 / 6000 [100%] (Sampling)
## Chain 2 finished in 94.5 seconds.
##
## Both chains finished successfully.
## Mean chain execution time: 93.4 seconds.
## Total execution time: 94.6 seconds.## Warning: 5 of 6000 (0.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.## Shape estimate: 9.720263
## Scale estimate 76.30104
## BRMS SUCCESS!
## BRMS RMSEs:
## $coef
## [1] 0.006483478
##
## $mrl
## [1] 5.998548
##
## $bioage
## [1] 7.433909
##
## **** COMPARISON ****
##
## **** SUMMARY COMPARISON ****
## Beta RMSE - AFT: 0
## Beta RMSE - BRMS: 0.0065
## MRL RMSE - AFT: 0
## MRL RMSE - BRMS: 5.9985
## BioAge RMSE - AFT 0
## BioAge RMSE - BRMS 7.4339## user system elapsed
## 316.119 27.228 116.8076.4 Method: Weibull; p = 200; \(n \in \{250, 500\}\); \(r_{nz} \in \{1\}\)
I will not do because takes too long!